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Abstract 

Current computational tools can generate and improve genome-scale models based on existing data; however, for 
many organisms, the data needed to test and refine such models are not available. To facilitate model 
development, we created the forced coupling algorithm, FOCAL, to identify genetic and environmental conditions 
such that a reaction becomes essential for an experimentally measurable phenotype. This reaction's conditional 
essentiality can then be tested experimentally to evaluate whether network connections occur or to create strains 
with desirable phenotypes. FOCAL allows network connections to be queried, which improves our understanding 
of metabolism and accuracy of developed models. 



Background 

There are currently over 3,000 completely sequenced 
bacterial genomes [1], For many of these sequenced 
organisms we know relatively little about them compared 
to well-studied organisms [2], even though they are 
important for biomedical, environmental, and biotechno- 
logical applications. However, their sequenced genomes 
provide a wealth of data that can be mined to discover 
their metabolic capabilities and transcriptional regulatory 
control mechanisms. Knowing how an organism metabo- 
lizes compounds, generates energy, produces cellular 
components, and synthesizes useful products is critical 
for enhancing chemical production, identifying new drug 
targets, or improving bioremediation. If little is known 
about an organism's metabolism and regulation a logical 
question is where to begin? Moreover, what sets of 
experiments should one perform to effectively determine 
how cells utilize and control metabolism? 

Mathematical representations of genome-scale net- 
works - known as genome-scale models (GEMs) - enable 
a quantitative and systematic approach to address this 
issue. By developing GEMs, the microbial reaction net- 
works can be interrogated to predict growth phenotypes, 
guide metabolic engineering strategies, elucidate network 
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components and interactions, and facilitate hypothesis- 
driven discovery [3-6]. However, the successful applica- 
tion of in silico metabolic and regulatory models depends 
on their ability to capture the underlying characteristics 
of the biochemical networks in the microbe of interest. 
With increasing improvements in genome sequencing 
technologies and annotation, and in metabolic network 
reconstruction [7], the ability to construct GEMs has 
become more high-throughput. Many of these annota- 
tion-derived GEMs possess reactions whose inclusion is 
based solely on homology or on reproducing growth phe- 
notypes (that is, enabling biomass production); conse- 
quently, verifying the metabolic networks derived from 
genomic data is becoming increasingly important. With- 
out an accurate representation of the microbial network, 
model driven design of therapeutics and metabolic engi- 
neering strategies will be potentially flawed and substan- 
tial time and resources may be wasted. Unfortunately, 
reactions and gene-protein-reaction (GPR) associations 
can be incorrectly included or omitted during model 
development due to database, sequencing, and annota- 
tion errors, as well as unknown enzyme functionality [4], 
Existing models for Escherichia coli have been painstak- 
ingly developed and refined over the past 20 years, using 
analysis of experimental data acquired over the past 50 
years from hundreds of laboratories. Spending this level 
of time, effort, and resources to obtain a good under- 
standing of metabolism for every microbial organism of 
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interest is simply intractable. Thus, to streamline the pro- 
cess of model curation, future experiments should be 
designed to reduce experimental efforts while still effec- 
tively probing the biological system of interest. Having 
the ability to quickly design experiments to test reactions 
is critical for improving the accuracy and utility of gen- 
ome-scale models, particularly for less-characterized 
microorganisms where existing experimental data are 
limiting. 

A GEM can be refined when discrepancies are found 
between model predictions and experimental observa- 
tions. Several automated computational approaches have 
been developed to suggest model improvements based 
on such discrepancies between model predictions and 
existing experimental data. Constraint-based model 
refinement algorithms, such as OMNI, SMILEY, Grow- 
Match, and GeneForce [8-11], work to improve a model's 
ability to reflect known experimental results. Depending 
on the algorithm, this may be accomplished by adding or 
removing network reactions, modifying GPR associations, 
modifying biomass compositions or relaxing regulatory 
rules. These methods successfully improve model accu- 
racy; however, they all rely on available experimental data 
to first identify model inaccuracies. 

Currently, there are no constraint-based methods to effi- 
ciently design new experiments to test the accuracy of a 
given genome-scale metabolic model, and its associated 
metabolic network reconstruction. To address this limita- 
tion, we sought to develop an approach that would identify 
media and gene knockout conditions under which a cho- 
sen reaction is essential for some measurable phenotype 
(for example, growth). A prior study has used minimal cut 
sets (MCSs) to identify minimal sets of reactions that if 
deleted will disable growth [12], and once enumerated 
MCSs could be evaluated to find a MCS involving the cho- 
sen reaction. However, identifying these sets requires com- 
putation of elementary modes, and so it can not be applied 
to genome-scale networks, which often contain approxi- 
mately 500 to 2,000 reactions [13]. Flux balance analysis 
(FBA) [6] can be used to predict if a reaction is essential 
for growth in genome-scale networks; however, finding 
conditions under which a chosen reaction is essential may 
require an exhaustive search of multiple gene knockout 
combinations. Additionally, since FBA predictions and 
MCSs are condition-specific, these methods would need 
to be evaluated in all possible media combinations, making 
the task even more computationally challenging. 

To address this experimental design challenge, we used 
concepts from flux coupling analysis to efficiently identify 
media and knockout conditions under which a chosen 
reaction is required to enable flux through another experi- 
mentally measurable reaction (for example, growth). Flux 
coupling analysis characterizes the relationships between 
reactions in a fixed network [14], and has been used to 



investigate gene regulation and gene essentiality [15,16], 
and for metabolic flux analysis [17,18]. In flux coupling 
analysis, all reversible reactions are first decoupled into a 
forward and reverse reaction. Then, the maximum and 
minimum flux ratio between two reactions is calculated 
and used to characterize the relationships between fluxes 
(v) through these two reactions. For example, if the mini- 
mum flux ratio (v c hosen/v m easured) is positive, then it implies 
that a chosen flux, v chosen , must be non-zero if another 
experimentally measurable flux, v measured , is non-zero 
(^measured — * v chosen)' For our purposes, reactions are con- 
sidered coupled if the minimum flux ratio is positive or 
the maximum flux ratio is a finite number; otherwise, they 
are uncoupled. These reaction couplings are highly depen- 
dent on the network and the environmental conditions 
used [14] and so flux coupling analysis has to be reapplied 
if the network changes (for example, a gene or reaction is 
deleted or added), or a different experimental condition is 
used (for example, glucose versus xylose media). As such, 
flux coupling analysis cannot identify network or environ- 
mental changes that lead to coupling between a chosen 
flux and an experimentally measurable flux. Thus, we 
developed the forced coupling algorithm (FOCAL) that 
will identify media conditions and gene deletions (which 
together form the coupling conditions) such that chosen 
fluxes are coupled with some measurable flux (that is, a 
flux that can be measured directly in experiments). Under 
these conditions, flux through a measurable reaction (for 
example, biomass production or by-product secretion) 
requires flux through one or multiple chosen reaction(s), 
and we refer to these conditions identified by FOCAL as 
coupling conditions. 

By finding coupling conditions in which biomass pro- 
duction depends on flux through a chosen reaction(s), we 
can design new growth phenotyping experiments to detect 
whether a chosen reaction occurs by simply monitoring 
cellular growth. Experimentally testing these coupling con- 
ditions allows for a variety of interesting conclusions to be 
made about the metabolic network. First, if no growth 
under the proposed coupling conditions occurs, then 
there is a problem with the model. In this case it is possi- 
ble that the chosen reaction does not occur because the 
associated enzyme is not expressed under this condition 
(due to regulation) or that the enzyme does not catalyze 
the reaction of interest (incorrect annotation). This means 
that regulatory, reaction and/or GPR changes are needed 
to correct the model. Second, if the chosen reaction is 
found to be conditionally essential under the coupling 
condition (meaning growth occurs under the coupling 
condition but when the chosen reaction is additionally 
eliminated no growth occurs), then the chosen reaction 
and its associated GPR relationships appear to be correct 
within the model. Third, if the chosen reaction is not con- 
ditionally essential under the coupling condition, then 
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components (for example, reactions or isozymes) are miss- 
ing from the network and can be suggested using compu- 
tational approaches [8,9]. 

A cycle of model testing and improvement can be estab- 
lished by iteratively using FOCAL to design experiments, 
conducting the FOCAL designed experiments, and adjust- 
ing the model when discrepancies between model predic- 
tions and experimental results are found (Figure 1). 

By enumerating and testing such coupling conditions, it 
is possible to not only confirm the presence of existing 
network components and interactions, but also to discover 
new interactions within the cellular network when the 
experimental results do not agree with model predictions. 
Additionally, since GEMS are powerful tools for enhan- 
cing biochemical production [19], we have also used 
FOCAL to design strains with complex and atypical phe- 
notypes, such as the concurrent utilization of multiple 
substrates by a single strain. By combining our novel 
experimental design algorithm with existing approaches 
for refining models [8-10], we envision an integrated com- 
putational and experimental platform (Figure 1) will be 
established that enables rapid development of highly accu- 
rate models and improved understanding of microbial 
metabolism across a wide variety of organisms, including 
those that are not well characterized experimentally. 




Figure 1 FOCAL refinement cycle. The model testing and refinement 
cycle is a three part process. First, FOCAL is used to design experiments 
where a particular reaction should be essential. The necessary mutants 
and media are prepared and growth phenotype experiments are 
performed. If any discrepancies are observed, the errors are corrected 
using various methods to suggest model improvements. These 
modifications can subsequently be tested further by designing new 
FOCAL designed experiments based on the refined model. 



Results and discussion 

FOCAL builds on concepts from the flux coupling frame- 
work [14], where the latter is capable of determining the 
relationships between two reaction fluxes given a fixed 
network and environment. Unlike the flux coupling frame- 
work, FOCAL actively works to create coupling within a 
network by selecting genetic and environmental condi- 
tions such that flux through a particular reaction (v c h OS en) 
becomes essential for another measurable flux (v measure d)- 
While a variety of different types of flux coupling relation- 
ships exist [14], FOCAL looks specifically for circum- 
stances under which the existence of a particular 
measurable flux, v measure d, implies the existence of flux 
through another reaction, v c h OS en (and, from contraposi- 
tion, no flux through v c h OS en implies no flux through 
Vmeasured)- Here, we discuss FOCAL's proposed solutions 
for coupling reactions to biomass in four genome-scale 
metabolic models. Using these results, we illustrate 
FOCAL's utility for systematically evaluating and refining 
metabolic models by comparing FOCAL predictions to 
new and existing experimental results. We further show 
how FOCAL led to the discovery of a new isozyme (YeiQ) 
for two reactions in glucuronate and galacturonate catabo- 
lism. Finally, we demonstrate the use of FOCAL to design 
more complex phenotypes, such as mutants that must 
concurrently utilize glucose and xylose in order to grow. 

Forced coupling algorithm: an illustrative example 

Using a small reaction network, we will first demonstrate 
how FOCAL works and how to interpret its results 
(Figure 2). FOCAL proposes minimal media components 
and knockout mutations (if needed) such that flux 
through the chosen reaction is required for biomass pro- 
duction. In the first example, FOCAL's objective is to 
design an experiment to test if the v 2 flux occurs. In the 
wild-type network (Figure 2a), biomass production (v bio ) 
and v 2 are uncoupled due to alternative ways of making 
the two biomass components, F and H (for example, 
using v 3 or, if metabolite G ex is in media, v 10 ). FOCAL 
indicates that coupling between v 2 and vt,i 0 (vbio — > v 2 ) 
can be obtained by using metabolite A as the sole mini- 
mal media component and deleting genes associated with 
v 8 (Figure 2b). FOCAL can also be extended to design 
substrate co-utilizing mutant strains as shown in Figure 
2c. To accomplish this, FOCAL looks for coupling condi- 
tions composed of minimal media specifications and 
gene deletions so that multiple reactions, in this case sub- 
strate transporters (vi and v 10 ), are required in order for 
the cell to grow (v bio — > Vi and v 10 )- To achieve this, 
FOCAL recommends deleting genes associated with v 3 
and V8 and using both metabolites A and G in the mini- 
mal media. The resulting mutant requires both metabo- 
lites A and G to produce biomass components F and H, 
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Figure 2 An illustrative example of FOCAL. FOCAL is first used to couple cellular growth (v blo ) with a chosen reaction flux (v 2 ). (a) In the uncoupled 
system, v 4 is coupled with v 2 (that is, v 4 > 0 implies v 2 > 0) but Vbio is not coupled with v 2 . (b) In the coupled case, Vbio is coupled with v 2 (Vbio -> v 2 ). 
Here, metabolite A ex is the only nutrient (no G ex ), and a gene associated with v 8 is deleted such that the upper pathway is required to synthesize 
metabolite F. Under these circumstances, flux through v bio requires flux through v 2 . Moreover, removal of v 2 (along with v a ) will result in a non-viable 
cellular mutant, (c) FOCAL can also be used to create substrate co-utilizing mutants where deletion of v 3 and v 8 requires the co-utilization of 
metabolites A and G in order to produce both biomass components, F and H. 



respectively. In some instances, alternative FOCAL solu- 
tions will exist and these can be found using additional 
integer cut constraints (see Materials and methods for 
details). 

Application to genome-scale metabolic networks 

To determine sets of experiments to test for all metabolic 
reactions in a network, FOCAL was applied to every 
reaction present in genome-scale metabolic networks for 
Escherichia coli [20,21], Bacillus subtilis [22] and Pseudo- 
monas putida [23] using biomass (that is, growth) as 
v measured- For each model, we specified sets of selectable 
carbon sources, nitrogen sources, electron acceptors, and 
additional nutrients that can be used to compose the 
minimal media (Additional file 1) and additional algo- 
rithm parameters (for example, maximum number of 
deletions; see Materials and methods for details). Based 
on FOCAL results, the reactions in these networks were 
categorized as coupled (a coupling condition could be 
found by FOCAL), uncoupled (no coupling condition 
could be found) or blocked (a reaction is incapable of 
carrying flux when all possible nutrients are provided) 
(Figure 3a). Each FOCAL proposed strategy was further 
evaluated based on the number of gene deletions 
required to achieve the desired reaction coupling 
between a metabolic reaction and biomass production 
(Figure 3b). Across the four models, a coupling condition 
was found for approximately 60 to approximately 85% of 
the unblocked reactions, and approximately 35 to 60% of 
these cases did not require any gene deletions, indicating 
that the media conditions alone were enough to couple 
the reaction to biomass (common deletions for each 
model can be found in Table SI in Additional file 2). For 
the iJR904 E. coli network, we also assessed how these 
reaction categorizations (that is, coupled, uncoupled, and 
blocked) were distributed across different metabolic sub- 
systems (Figure 3c) and how media components were 
used (Figure SI in Additional file 2). In E. coli, the cell 



envelope biosynthesis and the cofactor and prosthetic 
group biosynthesis subsystems contain a disproportionate 
number of blocked reactions. This is mainly due to the 
absence of many cofactors and prosthetic groups in the 
biomass reaction. Transporter, nucleotide salvage and oxi- 
dative phosphorylation reactions were the most difficult to 
find coupling conditions for, which may be attributable to 
redundant pathways, multi-functional enzymes, multiple 
isozymes or FOCAL simulation parameters. For E. coli, 
glucose, ammonia, and oxygen were the most frequently 
used carbon, nitrogen and electron acceptors utilized. 
Interestingly, the additional nutrients used in FOCAL 
designed experiments for E. coli and B. subtilis were quite 
different (Figure S2 in Additional file 2), likely due to 
differences in transporters between the two models. 
We also investigated if the gene deletions selected by 
FOCAL for the iJR904 E. coli network were close to the 
chosen reaction that becomes coupled with biomass. The 
shortest path distance between deleted reactions found by 
FOCAL and the chosen reaction was calculated for all pro- 
posed gene deletions associated with a single reaction (see 
Additional file 2 for details). For both a directed and 
undirected version of the metabolic network, the reactions 
that FOCAL deletes to achieve coupling tend to be closer 
on average (2.80 for the directed and 2.42 for the undir- 
ected network) than would be expected if deletions were 
selected randomly (4.99 and 3.90 for the directed and 
undirected network, respectively; in both cases P-value 
<le-10 using one-tailed f-test). 

Further analysis was done to investigate why FOCAL 
could not find coupling conditions for the 115 reactions in 
the iJR904 E. coli network that could not be coupled to 
biomass. These 115 reactions were subsequently re-evalu- 
ated with FOCAL using a higher gene deletion limit (up to 
10 gene deletions), adding more measurable reactions that 
FOCAL could use as v measurec | besides biomass production 
(by expanding the Coupling set, described in Materials and 
methods), and expanding the list of additional nutrients. 
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Figure 3 Various FOCAL statistics for genome-scale models, (a) Percentage of blocked, coupled, and uncoupled network reactions for each 
model evaluated, (b) Percentage of unblocked reactions from each model that require 0 to 5 deletions to become coupled with biomass. 
Reactions with zero gene deletions can be coupled solely by modifying the media composition. For all models, except U01366, the number of 
deletions is the number of necessary gene deletions. For U01366, additional isozyme deletions may be necessary (the total number of gene 
deletions needed for U01366 can be found in Figure S3 in Additional file 2). (c) Distribution of UR904 reactions belonging to a given coupling 
category (coupled, uncoupled or blocked) across metabolite subsystems. The percentage (left) or number (right) of reactions within a given 
coupling category that belong to a particular subsystem is shown. The fully coupled metabolic subsystem in (c) is composed of metabolic 
subsystems in which all reactions could be coupled to biomass, and contains the citric acid cycle, pentose phosphate cycle, nitrogen, pyruvate 
and methylglyoxal metabolism, purine and pyrimidine biosynthesis, anaplerotic, and putative reaction pathways. 



With these three changes, approximately a third of the 
previously uncoupled reactions were coupled by FOCAL 
to a measurable flux (Table 1). The remaining reactions 
could not be coupled for a variety of reasons. Around 40% 
of uncoupled reactions were involved in highly robust and 
interconnected pathways where reactions are catalyzed by 
the same multifunctional enzyme. For example, six reac- 
tions in the nucleoside salvage pathway (NTPP1-3 and 
NTPP5-7) dephosphorylate nucleosides and are all cata- 
lyzed by MazG, making it difficult to find coupling condi- 
tions that force one reaction to be essential while 
producing a viable mutant. Additionally, some reactions 
(approximately 10% of uncoupled reactions), based on a 
directed shortest path analysis, were not connected to bio- 
mass. For the remaining reactions, coupling conditions do 
not exist because they are involved in recycling metabo- 
lites, only participate in futile cycles, or have alternative 



Table 1 Comparison of FOCAL results for UR904 and 
mutant aerobic growth phenotypes 



Category 



Frequency 



Percentage of 
uncoupled 



Can be coupled using: 






More deletions 


8 


7.0 


More measurable fluxes 


8 


7.0 


More deletions/ measurable 


8 


7.0 


fluxes 






Expanded additional nutrient 


M 


12.2 


set 






Still cannot be coupled 






because: 






Highly robust/connected 


48 


41.7 


No connection to biomass 


11 


9.6 


Other reasons 


18 


15.7 


Total 


115 


100 
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reactions that cannot be eliminated due to their GPR rela- 
tionships (see Figure S4 in Additional file 2 for examples). 

Compared to the smaller E. coli model (iJR904), cou- 
pling conditions were found for a lower percentage 
(approximately 60%) of the unblocked reactions in the 
most recent E. coli model (iJ01366) [21]. We further 
investigated why coupling conditions could not be found 
for a larger fraction of these iJ01366 reactions, many of 
which involved transporters and membrane lipid metabo- 
lism (40% of all uncoupled reactions; Figure S3 in Addi- 
tional file 2). In many cases, no coupling conditions exist 
due to the presence of alternative reactions that are not 
associated with genes (for example, transporters like 
XANt2pp and XANtpp) or are associated with the same 
genes or essential genes. Each of the 24 EAR reactions, for 
example, has an alternative reaction that uses a different 
cofactor (NADPH versus NADH) and is associated with 
the same protein (FabI). As a result, the alternative reac- 
tions cannot be eliminated without also eliminating the 
chosen reaction. Other reactions that recycle metabolites 
back to their precursors were also in the uncoupled cate- 
gory since the recycling is never essential. The 98 phos- 
pholipase and lysophospholipase reactions that degrade 
phospholipids are examples of these. Another related pro- 
blem involves the irreversible export of compounds from 
the cytosol, which prevents their incorporation into bio- 
mass (for example, ZN2t3pp and ZN2abcpp), while other 
reactions cannot be coupled to biomass without adding 
compounds to the biomass equation. For example, the 14 
PSD and PSSA reactions produce phospholipids that are 
not part of biomass. 

Thus, an increased number of alternative reactions, recy- 
cling reactions and multifunctional enzymes in iJ01366 
reduces the number of reactions that can be coupled to 
biomass. As such, the increase in uncoupled reactions is 
not a failing of FOCAL, but rather a feature of the more 
comprehensive network. Future research could look to 
overcome this by instead generating coupling conditions 
for genes rather than reactions; in this way conditionally 
essential genes could be identified that would indicate that 
some of these uncoupled reactions take place. Addition- 
ally, while manual efforts were used to identify why parti- 
cular reactions cannot be coupled to biomass, this process 
could be semi-automated, by identifying clusters of reac- 
tions that share common genes and by determining cycles 
in metabolism (see Additional file 2 for details). 

Comparison of FOCAL predictions to experimental results 

FOCAL coupling conditions for E. coli iJR904 reactions 
associated with a single gene and involving only media 
specifications (that is, without requiring any gene dele- 
tions) were compared to previous studies where E. coli sin- 
gle knockout strains were tested for aerobic growth in 
glucose [24] and glycerol [25] minimal medium (Table 2). 



Table 2 Categorization of initially uncoupled reactions in 
UR904 





Glucose 


Glycerol 


Coupled reactions using only media 


193 


98 


Reactions associated with single genes 


13-1 


98 


Confirmed conditionally essential genes 


108 


44 


No experimental data 


0 


54 


Percentage agreement 


81% 


1 00% 



These experimental results were used to verify the condi- 
tional essentiality of the 232 single-gene reactions FOCAL 
coupled to biomass under these same media conditions. 
For the glucose experiment, a mutant was considered not 
to grow if the optical density (OD) at 24 and 48 hours was 
less than 0.10. For the glycerol aerobic experiment, we 
used the same growth classification as reported previously 
[25]. Of the 232 single-gene reactions that are coupled 
with biomass under these two conditions, experimental 
data were only available for 178 of the related mutants, 
and of these, 152 (approximately 85%) were conditionally 
essential, meaning that mutants missing these chosen 
reactions were unable to grow specifically under the pro- 
posed FOCAL media condition (Table 1). Of the 26 
model-data discrepancies, 2 mutants (AaroD and AnadC) 
were shown to be unable to grow on glucose in other 
experiments [26] and another 2 mutants (AfolB and AfolP) 
were shown to have gene duplications [27], indicating 
these 4 cases are likely not discrepancies. The remaining 
22 genes that were not conditionally essential indicate that 
changes to the model are needed. Model changes based 
on these datasets have been suggested previously [25] and 
involve: (1) eliminating components from the biomass 
equation; (2) accounting for additional transporters; and 
(3) adding isozymes or alternative reactions. This analysis 
illustrates how FOCAL results can provide confidence in 
model content and can lead to suggestions for improving 
the models when FOCAL predictions do not match 
experimental results. 

By determining coupling conditions for reactions with 
unknown GPRs, it is also possible to use FOCAL results 
to design high-throughput screens to identify genes 
associated with these so-called orphan reactions. Of the 
39 orphan reactions in iJR904 that are not transporters, 
coupling conditions were found by FOCAL for 27 reac- 
tions (Table S2 in Additional file 2). These coupling 
conditions can potentially be used to screen knockout 
mutant libraries to find conditionally essential genes 
that would be candidate genes responsible for these 
orphan reactions. The NAD-dependent succinic semial- 
dehyde reaction, SSALx (Figure 4), was one such orphan 
reaction, whose associated gene {ynel) has now been 
identified [28]. FOCAL predicts that the SSALx reaction 
is required in a AgabD mutant for aerobic growth on 
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Figure 4 Application of FOCAL to couple SSALx with biomass production, (a) The putrescine degradation pathway FOCAL predicts that 
the NAD-dependent succinic semialdehyde dehydrogenase (SSALx) reaction is coupled with biomass by removing gabD (whose gene product 
catalyzes SSALy) and growing the mutant in minimal media with putrescine and NH 4 + . (b) The feasible region for the hgabD mutant under this 
condition excludes points that lie on the x-axis with the exception of the origin (right), while the wild-type can grow on glucose without flux 
through SSALx (left). Metabolite abbreviations not reported in the text: a-kg, a-ketoglutarate; glu-L, L-glutamate. 



putrescine (ptrc) as a carbon source or on 4-aminobu- 
tanoate (4abut) as either a carbon or a nitrogen source. 
Both ptrc and 4abut are ultimately broken down into 
succinic semialdehyde (sucsal), which must be subse- 
quently consumed by one of the two succinic semialde- 
hyde dehydrogenases. The AgabD mutation prevents the 
NADP-dependent SSALy reaction from occurring and 
leaves NAD-dependent SSALx reaction as the sole suc- 
cinic semialdehyde dehydrogenase. Consequently, the 
AgabD mutant in one of these ptrc/4abut media condi- 
tions must use SSALx to grow and the desired flux cou- 
pling is obtained. Experimentally the AynelAgabD 
double mutant cannot grow during growth on putres- 
cine as a carbon source [29], indicating that FOCAL 
designed experiments can potentially be used to find 
genes for orphan reactions. 

To further illustrate the use of FOCAL in a model 
refinement cycle, growth phenotype experiments were 
performed based on FOCAL results for reactions in alter- 
native carbon metabolism. Reactions involved in galactur- 
onate and glucuronate catabolism (Figure 5) were selected 
due to the number of experiments proposed by FOCAL 
using these carbon sources and because reactions in these 
pathways were coupled to biomass using only media con- 
ditions allowing for facile testing (Table 3). All FOCAL 
predictions were consistent with measured single knock- 
out mutant growth phenotypes (that is, the genes asso- 
ciated with these chosen reactions were conditionally 
essential as predicted by FOCAL) with the exception of 



the AuxaB and AuxuB mutants, which grew on galacturo- 
nate and glucuronate, respectively (Table 4). Since the 
UxaB and UxuB enzymes carry out similar transforma- 
tions, we initially hypothesized that the two proteins may 
be able to catalyze both transformations. However, a dou- 
ble knockout AuxaBAuxuB mutant was still able to grow 
on both carbon sources. A BLASTp search found an oxi- 
doreductase gene with uncharacterized function, yeiQ, had 
significant homology to uxaB (E-value = e-21) and uxuB 
(E-value = e-155). Subsequent removal of yeiQ, uxaB, and 
uxuB eliminated the ability of strains to grow on glucuro- 
nate and galacturonate (Table 4). The results of these 
additional mutant phenotyping experiments (Table 4; Fig- 
ure S5 in Additional file 2) suggest that the altronate oxi- 
doreductase reaction could be catalyzed by UxaB, UxuB, 
or YeiQ and the mannonate oxidoreductase reaction 
could be catalyzed by UxuB or YeiQ. 

Substrate co-utilization strain designs 

FOCAL can also create more complex coupling condi- 
tions, where not just one but multiple reactions are 
coupled to biomass production. One such potential appli- 
cation of this is to design strains that co-utilize multiple 
substrates in order to overcome difficulties associated with 
diauxic growth and to speed up fermentation. Using this 
approach, a strain was proposed using the iJ01366 model 
for E. coli that is incapable of growth unless the cell con- 
currently consumes both glucose and xylose (Figure 6). 
This mutant has defects in both the pentose phosphate 
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D-galacturonate (galur) 

D-galacturonate isomerase 
{uxaQ 

D-tagaturonate (tagur) 
NADH 



altronate oxidoreductase 
{uxaB or uxuB oryeiQ) 



D-altronate (altrn) 

D-altronate dehydratase 
{uxaA) 




HX> 



D-glucuronate (glcur) 

D-glucuronate isomerase 
{uxaQ 

D-fructuronate (fruur) 
NADH 



D-mannonate oxidoreductase 
(uxuB oryeiQ) 



D-mannoate (mana) 

D-mannonate dehydratase 
{uxuA) 





2-dehydro-3-deoxy-D-gluconate (2ddglcn) 

Figure 5 D-galacturonate and D-glucuronate degradation pathways. Reactions involved in the degradation of galacturonate and 
glucuronate. Items in parentheses next to metabolites indicate metabolite abbreviations, and items in parentheses under enzymes indicate the 
associated genes. Gene names in black are those in the original UR904 model, while those in red indicate additional functionality discovered by 
FOCAL designed experiments that are added to the model to recapitulate experimental results. 



and glycolysis pathways, making it incapable of producing 
NAD/NADP and membrane lipids unless both glucose 
and xylose are consumed (see Table S3 in Additional file 2 
for a list of biomass components that cannot be made 
from individual sugars). The uptake of xylose and glucose 
concurrently allows the cell to produce dihydroxyacetone 
phosphate and glycerol-3-phosphate, which are used to 
produce NAD(P) and phospholipids. Such a mutant could 
be adaptively evolved to efficiently co-utilize both glucose 
and xylose under anaerobic conditions. 

A major distinction between this particular FOCAL 
designed mutant and others designed using elementary 
modes [30] is that we can consider genome-scale net- 
works and can enforce stricter co-utilization require- 
ments. Unlike previous designs that can utilize either 
glucose or xylose for growth and ethanol production 
[30], our algorithm identified a mutant where it is man- 
datory to use both glucose and xylose in order to grow, 
creating a strong selection for co-utilization in adaptive 
evolutionary experiments. Evolved mutants could 
improve lignocellulose conversion and avoid the added 
complications of developing and maintaining a co-cul- 
ture system [31]. By evolving co-utilizing mutants, pro- 
gress could be made towards more efficient strains for 
production of biofuels from lignocellulosic biomass. 



Conclusions 

FOCAL is capable of proposing experimental condi- 
tions (mutants and media composition) that will force 
coupling between a chosen flux of interest and a mea- 
surable flux (for example, cellular growth). As a result, 
FOCAL can design experiments to assess the accuracy 
and usage of metabolic reactions and their associated 
genes. FOCAL has numerous applications, including 
validating network elements, discovering new GPR 
associations and designing strains with unique and 
complex phenotypes. In addition, FOCAL coupling 
conditions could be used to select for improved 
enzyme activities since selection for improved growth 
would require improved flux through these reactions. 
Future work will look to reduce the total number of 
experiments needed to probe entire networks (by con- 
sidering alternative solutions) and incorporate more 
advanced modeling components such as regulatory 
information to improve strategies proposed by the 
forced coupling algorithm. 

Materials and methods 

Forced coupling algorithm 

FOCAL is a mixed-integer linear program (MILP) that 
works to propose media conditions and gene deletions 



Table 3 FOCAL designed experiments for reactions in galacturonate and glucuronate catabolism 



Chosen reaction 3 


Enzyme 


Associated gene 


FOCAL selected carbon source 


galur i± tagur 


D-Galacturonate isomerase 


uxaC 


D-Galacturonate 


h + nadh + tagur ^± altrn + nad 


Altronate oxidoreductase 


uxaB 


D-Galacturonate 


altrn -» 2ddglcn + h 2 o 


Altronate hydrolase 


uxaA 


D-Galacturonate 


glcur ?i fruur 


D-Glucuronate isomerase 


uxaC 


D-Glucuronate 


fruur + h + nadh ^ mana + nad 


D-Mannonate oxidoreductase 


uxuB 


D-Glucuronate 


mana -» 2ddglcn + h 2 o 


D-Mannonate hydrolyase 


uxuA 


D-Glucuronate 



Abbreviations match those shown in Figure 5. 
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Table 4 FOCAL experimental results 



Galacturonate 



Glucuronate 



Strain 



Experimental Model Experimental Model 



BW25113 

AuxaA 

AjxaB 

AuxaC 

AuxuA 

AuxuB 

AuxaB AuxuB 
AyeiQ 

AuxaB AyeiQ 
AuxuB AyeiQ 
AuxaB AuxuB AyeiQ 



NA 
NA 
+ a 
+ 

NA 



+ 

NA 

NA 



+ 
+ 
+ 

NA 



a Growth was delayed by >48 hours. NA, these experiments were not 
performed. 



such that a chosen flux (for example, fumarase) becomes 
coupled with another measurable flux (for example, cel- 
lular growth), meaning that flux through the measured 
reaction requires flux through the chosen reaction. 
FOCAL (summarized in Figure 7) is a bi-level algorithm 
composed of an inner problem that forces a flux ratio of 
interest to take its minimum value subject to media 
changes and gene deletions enforced by the outer pro- 
blem. The outer problem searches for media conditions 
and a minimal number of deletions such that the mini- 
mum flux ratio is positive, ensuring that coupling occurs 
between a measurable flux (from a user specified set, Cou- 
pling) and the chosen reaction (that is, v measure d — > v c h ose J- 
Since FOCAL calculates non-linear flux ratios, the pro- 
blem must first be linearized in order to solve the problem 
as a MILP. Therefore, the non-linear problem is trans- 
formed to its linear form as described previously [14], in 
this case by normalizing flux through each reaction (j), 
including the chosen flux, by the measured flux: 



^chosen 
^measured 



k 1 chosen -t — V chosen 



^measured 



(1) 



(2) 



For this transformation to be valid, all fluxes must be 
non-negative; thus, each reversible reaction was decom- 
posed into a forward and reverse reaction, and the resul- 
tant fluxes transformed as above: 



v i = v, 



%for 



Vj e R 



v ) t rev 



0, Vj ^ Irreversible 



(3) 



(4) 



where R is the set of all reactions, and R reV ersibie refers 
to the subset of all reversible reactions. FOCAL is for- 
mulated using the Equations 5 to 22 listed below. 



Outer objective 

max(r 0 f,j(A m£U . + 1) - 1) 

-aE x (l - ko g ) - PT,jmj,additional 

Inner objective 

min V c hosen,for + Vchosenjev 

Steady-state material balance 

E j € SijiVjjor - Vjjev) = 0, Vi € M 
[Unblocked) 



Uptake constraints 



MaxUptake 



, Vj e Exch 



hfor, Vj,™» f ^ 0 

Vj,rev — 0/ Vj ^ ^reversible 

Select measured flux 

Vj,f or = 1/ if ' Cj.for = L Vj e Coupling 

Vj.rev = 1/ ifcj.rev = L Vj € Coupling 

^ f ^ ^hfor + Cj,rev = 1 
Coupling 

Reaction deletions 

Vj,for,Vj,rev = 0, if Uj = 0, Vj £ R 

Media specifications 

v x ,rev = 0, if h x = 0, Vx e Exch\Minimal 

Coupling acceptance criterion 

{} 'chosen, for + Vchosen,rev) + (1 ^) — ^obj 

Media overlap rules 

m Xj k ■ media Xi it > h x , Vx e Exch 

keK 

ftix.k < hx, Vfe e K, Vx e Exch 
^keKmx.k ■ medial Xi k < 1, Vx e Exch 

Media uptake rules 

E x e m Xi k ■ media Xi k < maxMediai,, Vfe e k 

Exch 

GPR rules 

a i=f ( ko s) 



Deletion constraints 



(5) 

(6) 
(7) 

(8) 
(9) 
(10) 

(11) 
(12) 

(13) 
(14) 
(15) 
(16) 

(17) 

(18) 
(19) 

(20) 
(21) 
(22) 
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Figure 6 FOCAL predicted glucose/xylose co-utilization conditions, (a) Maximum FBA predicted anaerobic growth of the FOCAL designed 
£. coli strain as a function of the xylose fraction of the carbon source. The ratio of glucose and xylose within the minimal media was varied 
while maintaining a constant carbon uptake into the iJOl 366 network (1 10 mmol GgDW~ 1 -h~ 1 ). Under FOCAL's proposed conditions, the strain is 
incapable of growth when the media is composed entirely of glucose or xylose due to an inability to produce all biomass components. For 
comparison, the maximum predicted wild-type rate growth is 0.423 h" 1 on pure glucose and 0.362 h' 1 for pure xylose (not shown), (b) Possible 
fluxes through central metabolism in the mutant when grown only on glucose. Under these circumstances, the mutant is unable to produce 
dihydroxyacetone phosphate and glycerol-3-phosphate, which are critical for synthesizing NAD(P) and phospholipids. On xylose only (not 
shown), the mutant is incapable of sustaining flux beyond the pentose phosphate pathway. Boxed metabolites indicate biomass precursors and 
dashed arrows indicate multiple reaction steps. Metabolite abbreviations used but not provided in the text are: glc-D, D-glucose; g6p, glucose-6- 
phosphate; 6pgc, 6-phospho-gluconate; ru5p-D, D-ribulose 5-phosphate; r5p, ribose-5-phosphate; e4p, erythrose-4-phosphate; ftjp, fructoses- 
phosphate; fru, D-fructose; g3p, glyceraldehyde 3-phosphate; 13dpg, 3-phospho-D-glycerol phosphate; 3pg, 3-phospho-glycerate; 2pg, 2- 
phospho-glycerate; pep, phosphoenolpyruvate; pyr, pyruvate; accoa, acetyl-CoA. 



Inner primal problem 

The inner problem (Equations 6 to 15) minimizes the 
ratio of the two fluxes for both chosen reaction direc- 
tions such that if no coupling exists the inner objective 
is zero. This effectively amounts to solving the flux cou- 
pling framework problem proposed by Burgard et al. 



[14] to determine flux coupling, except, for their pur- 
poses, Burgard et al. also considered maximizing this 
objective. The transformed fluxes in the inner problem 
are subject to standard steady-state mass balance con- 
straints (Equation 7), which ensure that there is no net 
production or consumption for the set of all 



max ^ ^ ratio (Coupling Boolean) - Penalties 


(eq. 


5) 




subject to: min^ jt Linearized Coupling Objective 


(eq. 


6) \ 




subject to : Steady-State Material Balance (eq. 


7) 1 


Can be 


Uptake Constraints 


(eq. 


8-10) 


repeated for 


Measured Flux Selection 


(eq. 


11-13) 


multiple rxn 


Reaction Deletions 


(eq. 


14) 


couplings 


Media Specifications 


(eq. 


15) J 




Coupling Acceptance Criterion 


(eq. 


16) 




Media Overlap Rules 


(eq. 


17-20) 




GPR Rules (ko s -> aj) 


(eq. 


21) 




Deletion Constraints 


(eq. 


22) 





Figure 7 Overview of the forced coupling algorithm (FOCAL) 
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metabolites, M. Equations 8 to 10 are identical to those 
reported previously [14] to constrain substrate uptake 
and the linearization variable, t. Here, Exch is the set of 
all exchange reactions, and Vj is the maximal sub- 

strate uptake flux for that exchange (see Additional file 1 
for values used). Equations 11 to 13 are used to select the 
measured flux that will be coupled with the chosen flux 
of interest. The binary indicator variables, Cjj or and c^ rev , 
are used to determine which flux, from a specified set of 
measurable fluxes (Coupling), v chosen is being coupled 
with. Any deleted reactions (as determined based on 
GPR rules, described below), indicated by binary variable, 
fly, have their flux set to 0 in both the forward and reverse 
directions using Equation 14. All conditional constraints 
(Equations 11, 12, 14 and 15), were implemented using 
GAMS (GAMS Development Corporation, Washington, 
DC, USA) indicator constraints. 

To allow changes in minimal media conditions, four 
sets of media components were defined, each set con- 
taining exchange reactions used to import metabolites 
as sources of carbon, nitrogen, electron acceptors, or 
additional nutrients. Each media component set was 
specific for individual models. If experimental informa- 
tion was not available, FBA [6] was used to predict 
whether the microbe could use metabolites as a carbon, 
nitrogen or electron acceptor source or as an additional 
media component (media sets defined in Additional file 
1). Note these media component sets are largely a book- 
keeping mechanism for the user; thus, while a compo- 
nent may be categorized as a particular nutrient source, 
this does not mean that the organism will use this meta- 
bolite strictly for this purpose (for example, putrescine 
may be selected as nitrogen source, but may also be 
used as a carbon source in the model). Equation 15 
allows FOCAL to define the minimal media to be tested 
whilst removing all unselected substrate exchanges. 
Here, Minimal is the set of exchange fluxes that are 
essential for cellular growth irrespective of the carbon, 
nitrogen or electron acceptor selected (for example, 
water, protons, essential salts, ions, phosphate, and sul- 
fur sources). 

Outer problem 

In the outer problem, a binary indicator variable, r ob p is 
used to determine whether the desired coupling criter- 
ion has been satisfied (that is, V c hosen,for + V c hosen,rev >€) 

using the acceptance criteria constraint (Equation 16). 
For the E. coli and B. subtillis models, e was set to 1CT 5 
while for P. putida this was increased to 10" due to 
scaling differences between the models. To allow 
FOCAL to design media conditions, the media selection 
rules (Equation 17 to 19) were implemented as part of 
FOCAL's outer problem in which w,^ is a binary indi- 
cator variable used to select a metabolite exchange, x, 



from one of the created media component sets, while 
media^ is a binary matrix indicating whether metabo- 
lite exchange, x, belongs to the media component type, 
k. K is the set of four media component types (carbon 
source, nitrogen source, electron acceptor, and addi- 
tional nutrients), and h x is a binary variable used to con- 
trol the media composition and uptake rates in the 
inner problem (Equation 15). An optional constraint 
(Equation 19) prevents a given metabolite exchange 
from being selected as more than one media component 
type. Equation 20 also limits the total number of meta- 
bolite exchanges that can be used for each media com- 
ponent type. The parameter maxMedia^ was set to one, 
except for the co-utilization case, where it was set to 
two to enable use of two carbon sources. 

FOCAL is subject to GPR deletion rules (Equation 21), 
which were implemented as described previously [32]. 
Such rules use a series of binary variables to map gene 
deletions (ko g = 0) to associated reaction deletions (aj = 0). 
Limits on the maximum and minimum number of gene 
deletions were also imposed considering the set of genes 
in the model, G (Equation 22), using parameters A max and 
A min . For these studies, A max and A min were normally set 
to 5 and 0, respectively. In FOCAL's outer objective func- 
tion (Equation 5), a and ft. are positive real numbers used 
to penalize the use of gene deletions and metabolites from 
the additional nutrient set. For this study, values of a = 1.0 
and P = 0.25 were used so that adding additional nutrients 
would be favored over creating extra deletions, which take 
more time experimentally. FOCAL is not very sensitive to 
the penalty values, so these values can be changed to mod- 
ify the type of proposed experiments as long as the maxi- 
mum possible combined penalties do not exceed the 
increase in the objective resulting from the desired cou- 
pling. To solve the bi-level problem using available MILP 
solvers, the inner problem is rewritten using duality such 
that the primal and its dual are solved simultaneously and 
their objectives set equal to one another. This guarantees 
that the inner problem objective is met prior to maximiz- 
ing the outer objective [33]. Complete formulation of 
FOCAL as a single-level MILP is provided in Figure S6 in 
Additional file 2. An implementation of FOCAL in GAMS 
(GAMS Development Corporation) for the example net- 
work shown in Figure 2 is provided in Additional file 3. 

Evaluation of different networks 

Using FOCAL, coupling conditions were proposed for 
reactions within the genome-scale models of E. coli 
[20,21], B. subtilis [22], and P. putida [23]. Given the 
increased size of the network and complexity of GPRs 
in iJ01366 (Table S4 in Additional file 2), we first 
reduced the number of gene deletion decision variables 
for this model by excluding subunits and isozymes as 
described by Hamilton and Reed [34]. We also replaced 
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the Nuo and Ndh reactions in iJ01366 with average 
reactions since the flux through these reactions was 
constrained to be equal [21]. Additionally, we removed 
the wild-type biomass from the network and based all 
coupling off of the core biomass equation. 

For simplicity, a reaction flux was considered coupled 
to a measurable flux if a media and gene knockout strat- 
egy could be generated for either its forward or reverse 
component (Equations 6 and 16). To improve run-time 
performance, the set of possible measurable fluxes (that 
is, those in Coupling) that a chosen reaction could be 
coupled with initially only contained the biomass flux. 
For E. coli reactions for which FOCAL could not initially 
find a coupling condition, FOCAL was re-run using an 
expanded Coupling set that included ethanol, formate, 
and succinate secretion in addition to the biomass, since 
these metabolites are common anaerobic by-products 
and can be easily measured. 

CPLEX can take a significant amount of time to find and 
prove that a solution is the global minimum. Since we 
were mainly interested in finding FOCAL solutions for all 
reactions in the genome-scale networks and not necessa- 
rily finding the global minimum, we limited the time 
FOCAL could spend searching for a better solution; how- 
ever, this is not required if one desires to obtain a global 
solution. Using a CPLEX option (tilim), the algorithm was 
allowed only 3 hours to find a solution for any given reac- 
tion coupling problem. To further reduce the time spent 
solving for an optimal solution, once a feasible solution to 
the coupling problem was discovered, the algorithm was 
only allowed an additional 10 minutes to search for a bet- 
ter solution using the GAMS BCH facility. To minimize 
the number of different minimal media conditions pro- 
posed and to prune simple coupling problems, a reduced 
set of metabolite exchange reactions composed of glucose, 
ammonium, and oxygen exchanges as well as the entire 
additional nutrient set was used for the initial 10 minutes 
of solution time. If no solution was found within this time 
period, then a more exhaustive search was performed 
using all elements within the various media component 
sets for the remainder of the allotted 3 hours. This amount 
of time is comparable to other bi-level MILP methods 
given the number of decision variables involved. Further 
improvements in run-time performance may be possible 
by constraining the dual variables [35] or eliminating gene 
deletion decision variables for reactions that are coupled 
to other reactions under all media conditions [14]. Media 
component sets for the different models and run-time sta- 
tistics are provided in Additional file 1 and Table S5 in 
Additional file 2 respectively. 

Discovery of alternative solutions 

FOCAL will initially only propose a single coupling con- 
dition that best maximizes the objective. Under certain 



circumstances, alternative solutions may exist and can 
be found by adding integer cut constraints that make 
prior FOCAL solutions infeasible: 

(1 - ko s ) + (hj) < \OUKO\ + \OldMedia\ - 1 ^3) 

SeOMKO jeOldMedia 

where OldKO is the set of genes deleted in a past 
solution and OldMedia is the set of media components 
proposed in that same solution. Such a cut prevents 
FOCAL from proposing a solution that is identical to or 
a superset of a previous solution. Additionally, one can 
omit the knockout or media component of the integer 
cut depending on the type of alternative solutions one is 
interested in obtaining. 

Strains 

E. coli strains from the Keio collection [24], specifically 
uxaAy.kan, uxaBy.kan, uxaC::kan, uxuA::kan, uxuB::kan, 
yeiQy.kan, and E. coli K-12 BW25113, were used in 
FOCAL designed experiments. Additionally, three double 
mutants (AuxuB-.-.kan AuxaB, AyeiQ-.-.kan AuxaB, and 
AyeiQvk&n AuxuB) and a triple mutant (AyeiQ::kan AuxuB 
AuxaB) were generated using sequential removal of the 
kan gene using FLP recombinase [36] and PI transduction 
[37] followed with selection for kanomycin resistance. 

Growth phenotype plate experiments 

All strains were grown in triplicate at 37°C in a Tecan 
Infinite 200 microplate reader (Tecan Group Ltd, Swit- 
zerland) using 96-well plates. OD measurements were 
taken at 600 nm every 15 minutes with linear shaking 
(830 seconds, 4.5 mm). Tecan OD measurements were 
converted to an equivalent OD 600 value in a Biomate 
spectrophotometer with a 1 cm path length (see [10] for 
conversion factors used). All strains were pre-cultured 
for approximately 24 hours in M9 medium supplemen- 
ted with 2 g/L glucose and subsequently washed twice 
with M9 minimal media containing no carbon source to 
remove any residual glucose. Cells were then resus- 
pended in different media - M9 + 2 g/L D-galacturonate 
or M9 + 2 g/L D-glucuronate - such that the starting 
OD 60 o measurement was approximately 0.05 and then 
grown in the Tecan plate reader. 

Additional material 



Additional file 1: List of the maximum uptake rates and media 
components for the three genome-scale models 

Additional file 2: Supplementary material, including additional 
algorithm details, supplementary tables, and supplementary 
figures 

Additional file 3: An implementation of FOCAL for the example 
network shown in Figure 2 This file can be run by GAMS (GAMS 
Development Corporation, Washington, DC), which can be freely 
downloaded at [38]. 
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